1 Carga de paquetes

for (package in c("tidyverse","fpp3", "GGally", "normtest")) {
    if (!require(package, character.only=T, quietly=T)) {
        install.packages(package)
        library(package, character.only=T)
    }
}
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✓ ggplot2 3.3.5     ✓ purrr   0.3.4
## ✓ tibble  3.1.5     ✓ dplyr   1.0.7
## ✓ tidyr   1.1.4     ✓ stringr 1.4.0
## ✓ readr   2.0.2     ✓ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
## ── Attaching packages ──────────────────────────────────────────── fpp3 0.4.0 ──
## ✓ lubridate   1.8.0     ✓ feasts      0.2.2
## ✓ tsibble     1.1.1     ✓ fable       0.3.1
## ✓ tsibbledata 0.4.0
## ── Conflicts ───────────────────────────────────────────────── fpp3_conflicts ──
## x lubridate::date()    masks base::date()
## x dplyr::filter()      masks stats::filter()
## x tsibble::intersect() masks base::intersect()
## x tsibble::interval()  masks lubridate::interval()
## x dplyr::lag()         masks stats::lag()
## x tsibble::setdiff()   masks base::setdiff()
## x tsibble::union()     masks base::union()
## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2

2 Funciones auxiliares

# Contraste para los coeficientes
my_t_test <- function (object, ...) 
{
  par <- rbind(t_stat=tidy(object)$statistic, p_value=tidy(object)$p.value)
  colnames(par) <- tidy(object)$term
  if (NCOL(par) > 0) {
    cat("\nt-test:\n")
    coef <- round(par, digits = 4)
    print.default(coef, print.gap = 2)
  }
}

# Gráfico de correlogramas de residuos
my_tsresiduals <- function (data, ...) {
  if (!fabletools::is_mable(data)) {
    abort("gg_tsresiduals() must be used with a mable containing only one model.")
  }
  data <- stats::residuals(data)
  if (n_keys(data) > 1) {
    abort("gg_tsresiduals() must be used with a mable containing only one model.")
  }
  gg_tsdisplay(data, !!sym(".resid"), plot_type = "partial",  
    ...)
}

# Validación cruzada anidada

nested_cv <- function(df, h, last_train, string_formula){
  
  nested_errors <- vector()  
  
for (i in seq(last_train, last(df$date), h)){
  train <- df %>% filter(date<=i)
  test <- df %>% filter(date>i)
  
fitted_model <- train %>%
  model(arima = ARIMA(as.formula(string_formula)))

h_forecast = min(dim(test)[1], h)

fc <- fitted_model %>%
  forecast(h=h_forecast)

test_err <- fc %>%
  accuracy(test) %>%
  select(MAPE)
nested_errors <- c(nested_errors, test_err$MAPE)
NewList <- list("errors"=nested_errors, "mean"=mean(nested_errors))
}
return(NewList)}

# Test de autocorrelacines de los residuos
autocorrelation_test_plot <- function(aug, dof = 4, m = 12,  h = 5, alpha = 0.05){
vec <- c()
num_lags = seq(1, h*m)

for (i in num_lags){
 vec <- c(vec,aug %>% features(.resid, ljung_box, lag=i, dof=dof) %>% .$lb_pvalue)
}

autocorr_pvalues_resid <- tibble(
  lag = num_lags, 
  p_value = vec,
  incorelated = p_value >= alpha
)

plot <- autocorr_pvalues_resid %>% 
  drop_na() %>% 
   ggplot(aes(lag, p_value, color = incorelated)) + 
  geom_point() + 
  geom_hline(aes(yintercept = alpha), linetype="dashed", color = "indianred2")


newList <- list("values" = autocorr_pvalues_resid, "plot" = plot)
return(newList)
}

3 División train y test

Antes de empezar a hacer el estudio hay que distinguir la parte de los datos que se utilizará para construir la fórmula (para modelizar) de aquella que se utilizará para validar los resultados y cuyos datos, no deben ser utilizados.

Convertimos la serie en un objeto tsibble, para trabajar de manera cómoda. Analizamos visualmente la serie y las marcas de tiempo para saber cómo dividir los datos.

TU <- read_csv('SerieTU_format.csv')
## Rows: 156 Columns: 2
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl  (1): value
## date (1): fecha
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
# arreglamos formato de fecha
TU$fecha = tsibble::yearmonth(TU$fecha)
TU = as_tsibble(TU, index = fecha)

TU %>%
  autoplot(value) +
    labs(title = "Transporte urbano") +
    xlab("Fecha") + ylab("Miles de viajeros") +
  geom_point(aes(y = value), size = 1.2, shape = 16, color = "black")

Dado que tenemos datos mensuales y visualmente se aprecia una estacionalidad anual, reservamos el último año para test y el resto para train.

TU_train <- TU %>% filter_index(. ~ "2012-12-01")
TU_test <- TU %>% filter_index("2013-01-01" ~ .)

4 Fase de identificación

Se comienza la fase de identificación donde se deciden las transformaciones a realizar y el propio ajuste de la serie. Aunque veremos que cuando se realiza la selección de los parámetros del modelo ARIMA, la iteración con la estimación y el contraste son constantes hasta converger a un modelo válido.

4.1 Estacionariedad

Antes de analizar la serie mediante test estadísticos, analicemos gráficamente más en detalle, para confirmar la estacionalidad detectada.

TU_train %>%
gg_season(value, period = "year",labels = "left")

Estacionalidad anual: la forma de las series anuales son iguales.

Descomposición STL para analizar la existencia de tendencia y estacionalidad.

TU_train %>%
   model(seats = feasts:::STL(value)) %>%
  components() %>%
  autoplot()

Se confirman visualmente la estacionalidad anual y una tendencia creciente y luego decreciente.

A continuación, estudiamos si la serie es estacionaria en media y varianza y actuamos en consecuencia.

Estacionariedad en varianza

Analicemos qué transformación Box-Cox puede aplicarse para estabilizar la varianza.

lambda <- TU_train %>%
features(value, features = guerrero) %>% pull(lambda_guerrero)
lambda
## [1] 1.876598
TU_train %>% autoplot(box_cox(value,lambda)) +
labs(y = "Box-Cox transformed TU")

Dado que la serie visualmente muestra una varianza más o menos constante, que el valor Box Cox es más próximo a 1 que a otros valores de transformaciones usuales (logaritmo, raíz cuadrada) y que no ha habido un cambio muy llamativo con la transformación, no realizaremos ningún cambio a priori. Si en el proceso de estimación no conseguimos un modelo adecuado, volveríamos a este punto para valorar de nuevo la transformación.

Estacionariedad en media

La serie muestra una tendencia creciente y luego decreciente, es decir, la media va variando a lo largo de la serie. Por tanto, es esperable que se necesite realizar al menos una diferencia para estabilizar la media. Formalmente, se comprueba mediante una prueba de raíces unitarias.

TU_train %>%
  features(value, unitroot_kpss)

El p-valor es menor que 0.05, lo que indica que la hipótesis nula es rechazada. Es decir, hay evidencias estadísticas para rechazar la hipótesis nula de que la serie sea estacionaria. Por tanto, habría que realizar una diferencia regular.

TU_train %>%
  features(difference(value, 1), unitroot_kpss)

Una vez realizada la diferencia regular, el p-valor del test indica que no hay evidencias estadísticas para rechazar la hipótesis de que la serie sea estacionaria en media.

Graficamos la serie diferenciada y ya observamos que el nivel de la serie no cambia:

TU_train %>% autoplot(difference(value, 1)) + labs(title = "Serie con una diferencia regular")

Debido a la estacionalidad de los datos, quizá deba realizarse una diferencia estacional. Para evaluarlo podemos recurrir a la función unitroot_nsdiffs(), esta evalúa el estadístico fuerza estacional \(F_S\) y sugiere una diferencia estacional si esta es mayor que 0.64.

TU_train %>%
  mutate(turnover = difference(value,1)) %>%
  features(turnover, unitroot_nsdiffs)

El método sugiere una diferencia estacional. En cualquier caso, para asegurar que la serie es estacionaria, en la fase de identificación del modelo debemos asegurarnos de que las raíces de la parte autoregresiva están fuera del círculo unitario (y, por tanto, 1/raíz dentro del círculo).

4.2 Determinación del modelo

En los pasos previos se ha sugerido una diferencia regular y una estacional. Empezaremos teniendo en cuenta con la diferencia regular. Es decir, nos dirigimos a identificar los órdenes \(p, q, P, Q\) de un SARIMA (p,d,q)x(P,D,Q)\(_{12}\) donde ya sabemos que \(d=1\). Vemos que los residuos contienen mucha información.

fit <- TU_train %>%
  model(arima = ARIMA(value ~  1+pdq(0,1,0) + PDQ(0,0,0, period = 12)))
fit %>% my_tsresiduals(lag_max = 36)

En la ACF se aprecia una estructura autoregresiva regular y estacional cuyos órdenes se confirman en la PACF. Comenzamos introduciendo \(P=1\).

5 Fase de estimación y contraste

Comencemos ajustando la parte estacional mediante un SARIMA (0,1,0)x(1,0,0)\(_{12}\).

fit2 <- TU_train %>%
  model(arima = ARIMA(value ~ 1+pdq(0,1,0) + PDQ(1,0,0, period=12)))
fit2 %>% my_tsresiduals(lag_max =36)

report(fit2)
## Series: value 
## Model: ARIMA(0,1,0)(1,0,0)[12] w/ drift 
## 
## Coefficients:
##         sar1  constant
##       0.8615  -95.9125
## s.e.  0.0367  942.5549
## 
## sigma^2 estimated as 263289328:  log likelihood=-1596.33
## AIC=3198.66   AICc=3198.83   BIC=3207.55
my_t_test(fit2)
## 
## t-test:
##            sar1  constant
## t_stat   23.502   -0.1018
## p_value   0.000    0.9191

Hemos comprobado que el parámetro de este modelo es significativo y no existen problemas de estacionariedad.

Introducimos ahora \(q=1\). Es decir, SARIMA (0,1,1)x(1,0,0)\(_{12}\).

fit3 <- TU_train %>%
  model(arima = ARIMA(value ~ 1+pdq(0,1,1) + PDQ(1,0,0, period=12)))
fit3 %>% my_tsresiduals(lag_max =36)

report(fit3)
## Series: value 
## Model: ARIMA(0,1,1)(1,0,0)[12] w/ drift 
## 
## Coefficients:
##           ma1   sar1  constant
##       -0.8816  0.920  -21.7567
## s.e.   0.0428  0.025   64.3210
## 
## sigma^2 estimated as 118838387:  log likelihood=-1542.58
## AIC=3093.15   AICc=3093.44   BIC=3105
my_t_test(fit3)
## 
## t-test:
##              ma1     sar1  constant
## t_stat   -20.578  36.8724   -0.3383
## p_value    0.000   0.0000    0.7357

Introducimos \(Q=1\). Es decir, SARIMA (0,1,1)x(1,0,1)\(_{12}\).

fit4 <- TU_train %>%
  model(arima = ARIMA(value ~ 1+pdq(0,1,1) + PDQ(1,0,1, period=12)))
fit4 %>% my_tsresiduals(lag_max =36)

report(fit4)
## Series: value 
## Model: ARIMA(0,1,1)(1,0,1)[12] w/ drift 
## 
## Coefficients:
##           ma1    sar1     sma1  constant
##       -0.7892  0.9990  -0.8655    0.0141
## s.e.   0.0435  0.0017   0.1157    1.5368
## 
## sigma^2 estimated as 75708833:  log likelihood=-1520.22
## AIC=3050.44   AICc=3050.88   BIC=3065.26
my_t_test(fit4)
## 
## t-test:
##               ma1      sar1     sma1  constant
## t_stat   -18.1272  574.3631  -7.4833    0.0091
## p_value    0.0000    0.0000   0.0000    0.9927

El coeficiente estimado para “sar1” se acerca peligrosamente a 1. Es necesario realizar una diferencia estacional.

Diferencia \(D=1\). Es decir, SARIMA (0,1,1)x(0,1,1)\(_{12}\).

fit5 <- TU_train %>%
  model(arima = ARIMA(value ~ 1+pdq(0,1,1) + PDQ(0,1,1, period=12)))
fit5 %>% my_tsresiduals(lag_max =36)

report(fit5)
## Series: value 
## Model: ARIMA(0,1,1)(0,1,1)[12] w/ poly 
## 
## Coefficients:
##           ma1     sma1   constant
##       -0.8385  -0.9923  -115.6608
## s.e.   0.0448   0.1955    34.8462
## 
## sigma^2 estimated as 64684904:  log likelihood=-1377.47
## AIC=2762.94   AICc=2763.26   BIC=2774.44
my_t_test(fit5)
## 
## t-test:
##               ma1     sma1  constant
## t_stat   -18.7347  -5.0763   -3.3192
## p_value    0.0000   0.0000    0.0012

Introducimos \(p=1\). Es decir, SARIMA (1,1,1)x(0,1,1)\(_{12}\).

fit6 <- TU_train %>%
  model(arima = ARIMA(value ~ 1+pdq(1,1,1) + PDQ(0,1,1, period=12)))
fit6 %>% my_tsresiduals(lag_max =36)

report(fit6)
## Series: value 
## Model: ARIMA(1,1,1)(0,1,1)[12] w/ poly 
## 
## Coefficients:
##           ar1      ma1     sma1   constant
##       -0.3627  -0.7216  -0.8861  -157.0759
## s.e.   0.0937   0.0657   0.1301    58.5297
## 
## sigma^2 estimated as 64256797:  log likelihood=-1371
## AIC=2751.99   AICc=2752.47   BIC=2766.37
my_t_test(fit6)
## 
## t-test:
##              ar1       ma1     sma1  constant
## t_stat   -3.8717  -10.9771  -6.8126   -2.6837
## p_value   0.0002    0.0000   0.0000    0.0082

Introducimos \(p=2\) y quitamos constante. Es decir, SARIMA (2,1,1)x(0,1,1)\(_{12}\).

fit7 <- TU_train %>%
  model(arima = ARIMA(value ~ pdq(2,1,1) + PDQ(0,1,1, period=12)))
fit7 %>% my_tsresiduals(lag_max =36)

report(fit7)
## Series: value 
## Model: ARIMA(2,1,1)(0,1,1)[12] 
## 
## Coefficients:
##           ar1      ar2      ma1     sma1
##       -0.7184  -0.3967  -0.3891  -0.7876
## s.e.   0.1308   0.1117   0.1412   0.0879
## 
## sigma^2 estimated as 64713127:  log likelihood=-1368.38
## AIC=2746.76   AICc=2747.24   BIC=2761.13
my_t_test(fit7)
## 
## t-test:
##              ar1      ar2      ma1     sma1
## t_stat   -5.4916  -3.5525  -2.7562  -8.9572
## p_value   0.0000   0.0005   0.0067   0.0000
gg_arma(fit7)

5.1 Diagnosis de residuos

Para completar la fase de contraste evaluamos los residuos. Veamos primero el histograma.

aug <-fit7 %>% augment()

# Histogram
aug %>%
  ggplot(aes(x = .resid)) +
  geom_histogram(bins = 50) +
  ggtitle("Histogram of residuals")

Test media = 0

Contraste t-student para contrastar si la media es 0.

# Student's t-Test for mean=0
t.test(aug$.resid)
## 
##  One Sample t-test
## 
## data:  aug$.resid
## t = -0.69101, df = 143, p-value = 0.4907
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -1682.5268   810.8847
## sample estimates:
## mean of x 
## -435.8211

Como p-value > 0.05 no existen evidencias significativas para rechazar la hipótesis nula de que la muestra tenga media 0.

Test autocorrelaciones

A continuación comprobamos que los residuos están incorrelados.

# Ljung-Box autocorrelation
aug %>% features(.resid, ljung_box, lag=12, dof=4)

Como el p-valor es >0.05, podemos concluir que los residuos están incorrelados.

Test homocedasticidad

Para contrastar la heterocedasticidad se puede utilizar una regresión media-dispersión. Calculando los grupos de manera anual, dado que esa es la estacionalidad de la serie.

log_log <- aug %>% as_tibble() %>% 
  group_by(year(fecha)) %>% 
  summarize(mean_resid = log(mean(.resid+1)), std_resid = log(sd(.resid+1))) 
  
  summary(lm(std_resid~mean_resid, log_log))
## 
## Call:
## lm(formula = std_resid ~ mean_resid, data = log_log)
## 
## Residuals:
##        1        2        3        5 
##  0.05104  0.10698 -0.51048  0.35246 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)    1.319      0.718   1.837   0.2076  
## mean_resid     0.986      0.106   9.304   0.0114 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4466 on 2 degrees of freedom
##   (8 observations deleted due to missingness)
## Multiple R-squared:  0.9774, Adjusted R-squared:  0.9661 
## F-statistic: 86.56 on 1 and 2 DF,  p-value: 0.01136

El p-valor de log(media) es menor que 0.05, se rechaza la hipótesis de que sean homocedásticos.

Test de normalidad

Para contrastar la normalidad realizamos el test Jarque-Bera para evaluar la normalidad.

# Jarque Bera test
jb.norm.test(na.omit(aug$.resid))
## 
##  Jarque-Bera test for normality
## 
## data:  na.omit(aug$.resid)
## JB = 32.549, p-value = 0.001

Como el p-valor es menor que el nivel de significancia 0.05, se rechaza la hipótesis nula de normalidad de los residuos.

6 Fase de predicción

Finalmente evaluamos la capacidad predictiva del modelo. Calculamos las predicciones en el conjunto de test, calculamos sus errores y los comparamos con los residuos.

# Residual accuracy
resids <- fit7 %>% 
  accuracy() %>% 
  select(-c(.model, .type, ME, MPE, ACF1, RMSSE)) %>% 
  mutate(Evaluation='Training') 

# Forecasting
fc <- fit7 %>%
  forecast(h=12) 

test_err <- fc %>% 
  accuracy(TU) %>% 
  select(-c(.model, .type, ME, MPE, ACF1, RMSSE)) %>% 
  mutate(Evaluation='Test')

# Show errors together
bind_rows(test_err, resids) %>% select(Evaluation, everything())

Además de evaluar los errores con estas medidas globales, podemos evaluar los errores cometidos año a año visualmente. Así, podemos detectar outliers o efectos de calendario no detectados. Pudiendo así corregirlos con modelos más complejos mediante variables exógenas.

aug %>%  ggplot() +
  geom_line(aes(x = fecha, y = .fitted), color="navy") +
  geom_line(aes(x = fecha, y = value), color="gray24") +
  # geom_line(data=demandaGas_forecast, aes(x = index, y = value), color="red4") +
  ggtitle("SARIMA train fitted values") +
  xlab('Dates') +
  ylab('Demand') + facet_wrap(vars(year(fecha)), scales = 'free')

También podemos dibujar las predicciones con los intervalos de confianza (aunque en este caso no sean representativos, dada la no normalidad de los residuos). Evaluamos así también la capacidad de generalización del modelo, comparando las predicciones con el conjunto de test.

h <- dim(TU_test)[1]

demanda_plot <- TU %>% filter(fecha>last(TU_train$fecha)-14 & fecha< last(TU_train$fecha) + 12 )

fit7 %>%
  forecast(h=12) %>%
  autoplot(demanda_plot)

Como sabemos que los modelos SARIMA no predicen bien más alla del orden más alto disponible, en este caso 12, los errores en el conjunto de test más allá de 12 datos se dipararían.